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The physical properties of a classical many-particle system with interactions given by a repulsive 
Gaussian pair potential are extended to arbitrarily high Euclidean dimensions. The goals of this 
paper are to characterize the behavior of the pair correlation function g2 in various density regimes 
and to understand the phase properties of the Gaussian core model (GCM) as parametrized by 
dimension d. To this end, we explore the fluid (dilute and dense) and crystalline solid phases. For the 
dilute regime of the fluid phase, a cluster expansion of <?2 in reciprocal temperature (3 is presented, the 
coefficients of which may be evaluated analytically due to the nature of the Gaussian potential. We 
present preliminary results concerning the convergence properties of this expansion. The analytical 
cluster expansion is related to numerical approximations for g-z in the dense fluid regime by utilizing 
hypernetted chain, Percus-Yevick, and mean-field closures to the Ornstein-Zernike equation. Based 
on the results of these comparisons, we provide evidence in support of a decorrelation principle for 
the GCM in high Euclidean dimensions. In the solid phase, we consider the behavior of the freezing 
temperature Tf (p) in the limit p — > +oo and show Tj (p) — » in this limit for any d via a collective 
coordinate argument. Duality relations with respect to the energies of a lattice and its dual are then 
discussed, and these relations aid in the Maxwell double-tangent construction of phase coexistence 
regions between dual lattices based on lattice summation energies. The results from this analysis 
are used to draw conclusions about the ground-state structures of the GCM for a given dimension. 



I. INTRODUCTION 



The statistical mechanics of many-particles systems play a central role in the design and control of materials. The 
thermodynamic and transport properties of a material system are determined by the interactions among constituent 
particles, and these interactions are intimately related to the short- and long-range order in the system. Although 
many-body interactions certainly influence the determination of physical properties in these systems, one may generally 
obtain an accurate approximation to the underlying physics by considering pair interactions <p(r),r = ||x — x'||, 
between particles. Systems of particles intera cting via either inverse-power-laws (e.g., 4>{r) ~ l/r n , n € NJ^or hard 
core (<p(r) = +oo, < r < D; cf>(r) = 0, r > £)) 3 * 4 * 5 I pair potentials have been extensively studied in this regard. 

For interacting polymers, however, each of the spatially-extended macromolccules has a significant number of 
degrees of freedom. Modeling the pairwise monomer interactions in this system can rapidly become computationally 
prohibitive. Fortunately, the problem may be simplified immensely by considering instead the interactions among 
the centers-of-mass of the polymers. This assumption is equivalent to enforcing an effective interaction among the 
macromolecules. However, since it is certainly possible for the centers-of-mass of any two polymers to overlap, this 
effective potential must contain the the essential property of being bounded. The Flory-Krigbaum pair potential 
0fk(Y), introduced in 1950, provides the following form for the effective interaction between the centers-of-mass of 
two polymer chains: 6 

where V^g and V^oiv denote the volumes of a monomer segment and a solvent molecule, respectively, N is the degree 
of polymerization, R g is the radius of gyration of the chains, x is a parameter that controls the solvent quality 
[0 < x < 1/2 denotes a good (i.e., conducive to repulsion) solvent and x > 1/2 a poor one (conducive to attraction)], 
and P — l/(fcsT) denotes the reciprocal temperature scaled by Boltzmann's constant ks- 
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The form of ([TJ suggests that we consider as a general form of the effective pair potential for this system: 



(r) = eexp 



a 



(2) 



which is the pair interaction for the so-called Gaussian core model (GCM), originally introduced by Stillinger .^1 Here, e 
and a determine the energy and length scales, respectively for the system. The physical properties of this model have 
been well-documented up to three dimensions; it is known that the system may undergo a fluid-solid phase transition 
for sufficiently low temperatures (fc^T/e ~ 0.01), and within the solid-phase region there exists a FCC-BCC (d = 3) 
transition as the system passes from low density to high densityP Furthermore, the GCM displays re-entrant melting, 
in which the melting temperature T m (p) as a function of density p (= N/V) approaches in the limit p — > +oo; in 
other words, a crystal in the GCM at positive temperature can always be made to melt by isothermal compression.^ 
Similarly, T m (p) — » as p — > 0; this behavior follows directly from the reduction of the GCM to a system of hard 
spheres in this limit 

Despite the extent of research currently being pursued in this field, little is known about the high-dimensional prop- 
erties of the GCM. This gap in knowledge is in spite of the recent interest in the physics of high-dimensional systems 
of particles; for examp le, Torquato and Stillinger have previously examined the question of packing hard spheres in 
high dimensions! 9 * 10 ! 1 -^ This problem has applications to abstract algebra, number theory, and communications theory, 
where the optimal method of sending digital signals over noisy channels corresponds to the densest sphere packing in a 
high-dimensional space p21 Besides providing an improvement on the Minkowski lower bound on the maximal packing 
density in d-dimensional Euclidean space M. d , Torquato and Stillinger were also able to provide evidence for a decor- 
relation principle of disordered packings in high dimensionsP This principle states that as the dimension d increases, 
all unconstrained correlations vanish, and any higher-order correlation functions g^ n '(xi, . . . ,x„) may be written in 
terms of the number density p and the radial distribution function g 2 (r) within some small error. For equilibrium 
systems, this simplification implies that g 2 (r) approaches its low-density limit as its dimensional asymptotic limit; 
i.e., the high-dimensional behavior of 52 (r) is similar to its low-density behavior for any finite d. 

With regard to classical fluids, Frisch and Percus^have examined the Mayer cluster expansions for a pair-interacting 
system in high dimensions and have shown that for repulsive interactions, the series are dominated by ring diagrams 
at each order in particle density p. Resummation of the series leads to an analytic extension in density from which 
the second virial truncation remains valid at densities higher than the density at which the series diverge. Doren 
and Herschbacri^ previously have developed a dimensionally-dependent perturbation theory for quantum mechanical 
systems from which they draw conclusions about the energy eigenvalues in "physical" dimensions from the information 
obtained for values of d where simplifications in the behavior of the systems may occur. It is thus clear that one 
may obtain keen insight into the physical nature of a many-body system from an exploration of its high-dimensional 
analogs. 

As a result, our focus in the present study is on the phase properties of the GCM in arbitrary Euclidean dimension 
d. We make the preliminary disclaimer that when we henceforth speak of arbitrary dimension d, we imply a Euclidean 
geometry (d € N). It is significant to note that the GCM is ideal for this analytical study since it has the property 
that 4> E L p (R d ) Vp G [l,oo]. Therefore, <fi is absolutely integrable, and the Fourier and inverse Fourier transforms 
are uniquely defined and constitute an isometryP^ We utilize the following definition of the Fourier transform (FT) 
of a function /(x): 

/(k)= / expH(k,x)]/(x)dx, (3) 

where / denotes the FT of / and (k, x) = X^=i ^% x i denotes the inner product of two (real- valued) ei-dimensional 
vectors. Similarly, the inverse FT is defined as: 

/(x) = (i) / Rd cx PK k . x )]/>) dk - ( 4 ) 

For radial functions [i.e., /(x) = /(||x||) = f(r)], fj) and @ take the forcnPl 



f(k) = Wl "^m ^^ dr (5) 
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The FT cf>(k) of the pair potential <p(r) in the GCM is given for all d by: 

$(k) = (7rcr 2 ) d/2 eexp 
and the integral of (f> over M. d is (for x,x' e R d ): 



/ 0(||x-x'||)dx = e(7TCT 2 ) d / 2 . (8) 

We immediately see that the dimensionality of the problem is contained entirely in the factors (na 2 ) d / 2 , facilitating 
the generalization of the model to arbitrary dimensionality Our goal is to characterize the fluid and solid phases of 
the GCM in high dimensions. 

With regard to the fluid phase, Stillinger and coworkers have developed high-temperature expansions for the excess 
free energy^ f(/3) and radial distribution function^ (72(f), the former for arbitrary d and the latter for d = 3. The 
convergence properties of the free energy expansion have been previously explored,^ although the series for /(/?) 
is divergent, it may be formally evaluated via use of Borel resummation. However, similar convergence properties 
for the (72(f) series remain unestablished even for d = 3. Furthermore, while compariso ns in three dimensions have 
been drawn between "exact" representations of g%{r), either from molecular simulations^^! r the aforementioned 
expansions,^ and numerical approximations, little is known of the high-dimensional applicability of numerical meth- 
ods. Of particular interest is the validity of the "mean-field appro ximati on" (MFA) to the direct correlation function 
c(r) that arises from a density functional description of the GCMPSEH By extending the temperature expansion of 
(72(f) and generalizing to arbitrary dimension, we attempt to elucidate the relationship among the MFA, hypernetted- 
chain (HNC) approximation, and the Percus-Yevick (PY) approximation to the GCM and relate the results to a 
decorrelation principle. Simultaneously, we explore the convergence properties of the high-temperature expansion of 
(72(f) in arbitrary dimension and relate the results to the phase behavior of the model. 

Prior work on the solid phase of the GCM involving lattice summation energies calculated in d — 3 for simple 
cubic (SC), body-centered cubic (BCC), face-centered cubic (FCC), hexagonal close-packed (HCP), and diamond 
(DIA) structures indicates a transition in the minimum energy of the system at p w 7r -3 / 2 from FCC to BCC.^ This 
conclusion has been supported and expanded to d — 1,2 via the calculation of duality relationships rela ting t he energy 
per particle (<j>/iV) A at low density of a lattice A to the corresponding ($/iV) A , of the dual lattice A* l 2 ° l 21 l However, 
there remains an open question of the relative stability of lattices in higher dimensions with respect to minimization 
of the lattice summation energy. This problem becomes especially apparent for d > 6, where the family of lattices 
to which FCC belongs no longer represents the densest known sphere packing among lattices.^ Worthy of mention 
in this regard is the corresponding conjecture by Torquato and Stillinger 22 that the Gaussian core potential and any 
other sufficiently well-behaved completely monotonic potential function share the same ground-state structures in M. d 
for 2 < d < 8 and d — 24 although not necessarily at the same densities; more specifically, they claim that these 
ground states are the Bravais lattices corresponding to the densest known sphere packings for < p < p\ and the 
corresponding reciprocal Bravais lattices for P2 < p < +00, where p\ and pi are the density limits for the phase 
coexistence region of the lattices!^ We seek to provide numerical support for the latter part of this conjecture with 
respect to the GCM. 

To these ends, we begin in Section II by developing the requisite high-temperature cluster expansion in (3 for 
ln[<72(r)] and explore the convergence properties of the series for arbitrary d. At low densities for which the series 
is appropriate, information about the behavior of the dilute fluid regime of the GCM may be thus obtained. In 
Section III we explore the dense fluid regime of the GCM using the three numerical approximations listed above in 
order to obtain information about (72(f) and the associated structure factor S(k) for the system. The validity of 
these approximations for arbitrary d is then evaluated. It may be shown that S(k) approaches a step function with 
discontinuity at k — +00 in the infinite-dimensional limit; we use this information to provide analytical support for 
a decorrelation principle in the fluid phase of the GCM. We devote Section IV to the solid phase of the GCM. The 
behavior of the melting temperature T m {p) in the limit p — *• +00 is generalized with respect to d, providing evidence 
for a fluid-solid phase transition at high density and sufficiently low temperature. We then calculate lattice summation 
energies for the A^, Dj, E^, and 1 d lattice families and their duals in various dimensions. After establishing duality 
relationships for these lattice families, we explore the phase coexistence regions between the lowest-energy lattices and 
their duals via Maxwell double-tangent constructions and show that the width of these regions increases with respect 
to the self-dual density p* as the dimensionality increases. The information gathered from this analysis provides 
evidence for the Torquato-Stillinger conjecture^! mentioned above concerning the ground states of certain classical 
many-particle systems. Concluding remarks are given in Section V. 
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II. DILUTE FLUID-PHASE VIRIAL BEHAVIOR OF THE GCM 

A. Correlation function formalism 

Although there is no evidence for a conventional gas-liquid phase transition in the GCM, it is mathematically con- 
venient to consider the "dilute" and "dense" fluid regimes separately; the reasons for this distinction will become clear 
momentarily Our understanding of the dilute (i.e., low density) fluid phase will involve the analytical determination 
of the radial distribution function g 2 (»") in terms of an infinite series in (3. To motivate this correlation function, we 
recall that for a system of N (fixed) particles, the configurational part of the canonical partition function is given by: 

Z N = ■■■ exp[- ( 3$(x 1 ,...,x w )]c2x 1 ---dxiv, (9) 
J Jv N 

where ^(xi, . . . , Xjv) = 4>{ r ij) denotes the total potential energy in the system. As a result, the probability 

l<i<j<N 

distribution for observing the many-body system with configuration x N = {xi, . . . ,xjv} is: 

P(x N ) = exp[-^(x N )]/Z N . (10) 



Based on (10 1, we define the n-body correlation function /?(™ ;JV )(xi, . . . ,x„) to be the joint probability distribution 



function that in an TV-particle system, particles will be found at positions {xi, . . . , x„}. Mathematically, 

^■ JV )(x 1 ,...,x„) = — - - [■■■ [ P(x N ) dx n+1 ... X2 v, (11) 

{N-n)\J J V (N-n) 

where the prefactor Nl/( N — n) ! denotes the number of ways of choosing an ordered subset of n particles from a total 
population of N. For an isotropic fluid (as found in the GCM), we have the result that: 

p^ N \ Xl ) = N/V = p. (12) 

It is therefore reasonable to introduce the functions g n (xi, . . . ,x„), defined by: 

/^(xx, . . . ,x„) = p(" ;Ar )(xi, • • . ,x„)//A (13) 

It is significant to note from the formulation in (13) that when j'"'(xi, . . . , x^) = 1, there is an absence of 



correlations in the system. Our primary concern in this study is with the function (/^(xi, X2), which for an isotropic 
fluid takes the form g 2 (r) with r = ||xi — X2 1| - The function 32 is known as the radial distribution function or pair 
correlation function; since 

— ~ = P92(r), 14 

P 

we identify 52 as being proportional to the conditional probability density that a particle will be found at radial 
distance r given that another is at the origin. Equivalently, pg2{f) is the average particle density at radial separation 
r given that a particle is located at the origin. 

Since for a general fluid the correlations in the system will diminish with increasing radial separation r, we have 
the asymptotic behavior g 2 (r) — > 1 as r — > +00. It is conventional also to introduce the so-called total correlation 
function h(r), defined by: 

h(r) = g 2 (r)-l. (15) 
It follows from the properties of g 2 that h(r) — > as r — > +00. 

B. Mayer cluster expansion of 

It has been well-documented that the radial distribution function g 2 may be expanded as an infinite series in the 
density p; this expansion has the formPSl 



g 2 {r) = exp[- 




En/- 



G m+2 \a£G„ 




(16) 
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where G m +2 denotes the graph containing m integrable vertices and two stationary vertices which becomes biconnected 
when an edge is added between the two stationary vertices. The parameter a = (i,j) corresponds to an edge in the 
graph G m+ 2, and f a — exp [— /3<j) a ] — 1 denotes the Mayer /-function with <p a representing the pair potential governing 
the interaction between particles located at Xj and Xj. A biconnected graph is a collection of vertices and corresponding 
edges such that one may trace a path between any two vertices even upon the removal of an edge. We define an 



integrable vertex to be any vertex in the graph G with a corresponding variable of integration Xj in ( 16 1; a stationary 



vertex has no corresponding variable of integration. Unless otherwise stated, vertices 1 and 2 will always be defined 
as the stationary vertices of any graph G m +2, meaning that §2 remains a function of r = ||xi — X2||. For example, the 
expansion of 172 up to 0{p 2 ) is: 



g 2 (r) = exp [-{3<p 



1,21 



1 



/l,3/2, 3 dx 3 + 0(p 2 ) 



(17) 



We note that in the limit p — > 0, g-z{r) — > exp [— (3<j>(r)\ , which is the Boltzmann factor for the pair potential <f>. For 
the GCM, it is more useful to pass to an expansion of gi in reciprocal temperature /3; this conversion is accomplished 
by a Taylor expansion of each of the /„ as follows: 



fa 



exp[ 

: (3(j) a 



(18) 
(19) 



Subsequent multiplication of the /-functions and recollection of orders of (3 leads to the desired expansion. The reason 
for passing to the reciprocal temperature expansion is that each of the integrals in ( 16 ) reduces to an integral over 



products of Gaussians, which may usually be evaluated analytically by repeated use of the relation (for e — a = 1): 



d/2 



,J i.k 



(20) 



Equation ( 20 1 follows directly from (|8| . 

Although, as mentioned above, it is possible to analytically evaluate each of the terms in the (3 series derived from 
(16 1, the mathematical complexity of the problem increases significantly around 0(/3 6 ) due to the increasing number 



of integrals to evaluate. However, it turns out that several of the graphs from this expansion may be "removed" by 
passing to the expansion of In [52 (V)] , which may be suitably derived from ( 16 1 by a Taylor expansion. The advantages 
of the logarithmic expansion are that (a) it contains the entire exp [— /3<^>x ,2 J term from (16 1 in —/3<pi t 2 and that (b) 



it removes all parallel graphs from the expansion, thereby drastically reducing the computational cost of analytical 
evaluation. In general, this series takes the form: 



+00 



In [ga(r)] = £ (-/?)"/»(»•)• 



(21) 



An investigation of convergence properties of the /3-series in (21 1 may be done via reference to the ratio test for 
infinite series We consider the ratios defined by: 



(n+i(r) = 



fn+i(r) 



Mr) 



fn+i(r) 
Mr) ' 



The radius of convergence A of the series is then given by: 



(22) 



limSUpCn+lM = y. 

n—*oo A 



(23) 



We note that the left-hand side of (23) is guaranteed to exist W e [0,oo) although it may be ±00. It is important 



to note that the proof of convergence of the density expansion (16 1 is well-known and can be found, for example, in 
RuellePS All that is required is to show that the GCM pair interaction is stable and regular. Borrowing the definitions 
from RuelleJ^a fc-body interaction (4>^)k>2 is stable if there exists a B > such that: 



for all n > and xi, 
and satisfies: 



$( Xl ,...,x n ) > -nB (24) 
, x„ 6 M. d . A pair interaction (f> is regular if it is bounded from below by a finite constant K 



C(f3) = / |ex P [-/fy(x)] - l|dx < +00 



(25) 
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for some /? > and hence for all (3 > 0. The proof that the GCM pair potential is regular can be found in Appendix 
A. 

We have successfully evaluated each of the f n (r) in (21 1 up to 0[(— /3) 5 ) with the results (in units such that 
e = o- = l): 



fi(r) = 



(26) 



Mr) = ,(ff> 



2\d/2 



(27) 
(28) 



2p 2 



d/2 



/,3/5 



d/2 



a/4 



(29) 



/ B (r) = 



2p 2 



3p J 



2\ <V 2 



3\ d/2 



k5/8 



3 / 8 + 2p J 



2p 3 



3\ <V 2 
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b 2 ' 7 + p 2 



1\ /7T 



2\ d/2 



(30) 



d/2 



a/ 5 



3 \ d / 2 



3\d/2 



a/ 2 



where, as above, 



= </>(ri2) denotes the pair interaction between two particles in the GCM. Using the results in 
(26l-(30l, we may obtain an approximate plot of 52 m the limit of high temperature and low density, indicative of a 
dilute fluid. Such a plot is included in Figure [T| for various values of d. 

Note that the value of <?2(0) is nonzero in each case; we expect this result since the bounded nature of the potential 
allows for a finite probability of particle overlap. Furthermore, the fact that the potential is repulsive suggests the 
presence of negative correlations near the origin, which we also observe. Correlations in the dilute fluid phase rapidly 
diminish as r increases; what is perhaps most significant about 32 in this case, however, is that the correlations also 
diminish as the dimensionality of the system increases; i.e., 32 ~~ * 1 more quickly as d increases. This observation 
provides the first suggestion (but of course does not prove) that a decorrelation principle applies for the GCM. We 
will return to this point in the study of the dense fluid regime. 



To help us understand the high-dimensional behavior of 52, we examine more closely the expansion in (21 1. Without 
loss of generality, we will momentarily work with reduced units such that a = e = 1. Let C n +\ denote the unique chain 
diagram contribution (containing n + 1 total vertices) to the order (3 n factor in the gi high-temperature expansion; 
namely, 



~ n n 
i=l i=2 



(31) 



where for notational convenience we have chosen vertices 1 and n+ 1 as the stationary vertices in the graph. We show 
in Appendix B that: 



C n +\ — p" 



1\ d/2 



a/ r 



(32) 
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FIG. 1: The radial distribution function for p = 0.1, (3 = 0.1 as obtained from the summation in (21 1 to 0[(— /3) 



Vn > 2. 



The argument used to prove (32 1 elucidates a central property of the cluster integrals: the exponential order of tt 



is determined solely by the number of integrable vertices in the corresponding cluster diagram. This statement is not 
exclusive to the chain diagrams. Since the dimensionality d is contained only in the factors (7r Q /^) d ^ 2 that appear 
for each cluster integral, we see that those diagrams that maximize the value of a for 7 of order n will dominate the 
<72 cluster expansion. Note that for n € N \ {1}, 7r n_1 > n, which may be easily proved by mathematical induction. 
Since the chain diagram for the (3 n contribution to the cluster expansion contains the greatest number of integrable 
vertices, we expect this diagram to maximize the quotient (n a /j) and thus dominate the expansion at any order 
(3 n for 7 of order n. This conclusion is in agreement with previously reported results for classical pair-interacting 
repulsive fluids by Frisch and PercusP^ Therefore, as the dimension d increases, we expect that 32 may be represented 
to a good approximation by the summation over these chain diagrams; we will see momentarily that this truncation 
of the series corresponds to the MFA and that the result is convergent. 

This convergence claim is not true, however, for the general series representation, as may be seen from the plots 
of the ratios ( n +i(r), defined in (22), in Figure [2] We note that the ratios depend strongly on the dimension d 
of the system and the radial separation r; however, they eventually diverge as r increases. There is no reason a 
priori why the gi cluster expansion should converge only for some values of r, and we therefore expect the series to 
diverge as a whole. We recall that the high-temperature expansion of 0.1 is centered at f3 = 0; therefore, if it is to 
converge for some j3 > 0, it must also converge for certain values of /3 < 0. However, as the reciprocal temperature 
(3 passes through 0, the pair potential undergoes an effective transition from repulsive to attractive interactions due 
to the coupling of (3 and the energy scale e in the Boltzmann factor. Because the pair potential does not exclude 
particle overlap, even in the limit r — > 0, the attractive regime forces the particles to cluster on top of each other at 
a point, and the system undergoes a collapse instability. This collapse instability forces the radius of con vergence to 
and has been documented by Stillinger for the equivalent high-temperature free-energy expansion! 17 ! 25 ! In terms of 
interacting polymers, the collapse instability corresponds to changing the effective composition of the solvent such 
that aggregation of the macromolecules is energetically favorable. The FK potential in ([TJ captures this behavior via 
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FIG. 2: The ratios C»i+i( r ) derived from (26l-(30l for p = pn d ^ 2 = 0.9. Reduced units e = a = 1 have been used 



variation in the parameter \- We make note, however, that the divergence of Cn+i( r ) appears to push outward toward 
r = +oo as d increases. It is therefore possible that a high-dimensional approximation such as the one mentioned 
above may be able to overcome the collapse instability, and we explore this possibility in Section III below. 



III. DENSE FLUID-PHASE BEHAVIOR OF THE GCM 

The low-density cluster expansion studied in the dilute fluid regime, though convergent in p according to Ruellej^ 
has a finite radius of convergence which is necessarily small. Therefore, to discern information about the fluid phase at 
values of p greater than the radius of convergence of the series above, we rely on approximation methods (numerical 
and analytical) to estimate g^- These approximation methods rely on solutions to the so-called Ornstein-Zernike 
equation, given by: 



h(r) = c(r) + p [ h(\\x - x'||)c(x')dx' 

JR d 

— c(r) + p(h * c)(r), 



(33) 
(34) 



where c(r) denotes the direct correlation function, and {h*c) indicates the convolution of h and c. We take the Ornstein- 
Zernike (OZ) equation as the definition of the direct correlation function; it essentially separates the immediate 
interactions between two particles from those interactions that result indirectly from interactions with surrounding 
particles. The advantage of the OZ equation is that it allows us to approximate h and therefore by making a 
reasonable ansatz about the form of c. We consider the following three well-documented closures to the OZ equation 
(see, e.g., McQuarrie^ 1 



CHNc( r ) = exp [— [3<j){r) + 7(r)] — "f(r) — 1, (hypernetted chain) 

cpy{t) = {1 + j(r)} {exp [— /3</>(r)] — 1} , (Percus-Yevick) 
cmfaM = — (3<fi(r), (mean-field approximation), 



(35) 
(36) 
(37) 



where 

~f(r) = h(r)-c(r). (38) 
Both the PY approximation and the MFA may be obtained from the HNC approximation via linearization of one or 



more of the exponential functions in ( 35 ) . Solutions to the OZ equation utilizing the HNC and PY approximations 
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are necessarily numerical. We utilize a relatively simple modified Picard iteration algorithm, a review of which may 
be found in the article by Busigin and Phillips!^ For reference, we include the details of the algorithm in Appendix 
C. 

The MFA is unique in the sense that one may obtain some analytical results from the OZ equation using Fourier 
analysis as a result of the relation Q. Namely, the OZ equation implies: 



h(k) = 

/iMFA(fc) = 



1 - pc(k) 



—f3e (7rcr 2 ) d ^ 2 exp 


(ak) 2 ' 
4 




1 + /3ep(7rcr 2 ) d / 2 exp 


( CT fc) 2 " 
4 



(39) 
(40) 



Our study of this approximation will involve the associated structure factor factor S(k) for the system at a given 
density, which we introduce here as: 



S(k) = l + ph{k) 



1 - pc{k) ' 



(41) 



where the second equality follows from the OZ equation. The structure factor is proportional to the scattered intensity 
of radiation from a system of points and thus is experimentally observable; this notion is the physical motivation for 
the intrinsic property that S(k) > Vfc G [0,+oo). Determination of the structure factor for a system is therefore a 
means to establish the physicality of a given approximation. For the MFA, note that (|7| and (37 1 imply: 



S'MFA(fc) = 



1 + (3ep(-K(j 2 ) d / 2 exp 



4 



(42) 



We collect in Figure [3] the results for gi as obtained from the HNC and PY approximations along with the MFA 
utilizing the iterative Fourier algorithm above. We have chosen the density p to be sufficiently high to capture the 
behavior of the dense fluid regime. Each of the approximations shows short-range correlations which diminish rapidly 
with increasing r, indicative of gas-like fluid behavior. Again, the probability of finding particles at zero separation 
is nonvanishing due to the bounded nature of the potential. However, while the HNC approximation and the MFA 
have similar values for 52(0), it appears that the PY approximation underestimates this value and thereby introduces 
an increased effective repulsion among the particles. Note that the correlations also diminish with dimensionality, 
providing further numerical support for a decorrelation principle with respect to gi- 

Figure [4] shows the results for the HNC approximation with d = 5 and (3e < 0. Here, the approximation does 
not immediately diverge and shows aggregation at the origin; as the value of (3e becomes increasingly negative, the 
correlations in the system become longer in range until the approximation does in fact diverge for j3e slightly less than 
—0.73. Similar results have been prior ly reported by Root, Stillinger, and Washington 18 for the PY approximation. 
Thus, these numerical schemes have the capacity to branch into the instability region for at least some small range 
of f3e < 0. The results obtained in this region certainly do not contain the actual physics of the GCM, but they do 
reflect the properties of aggregation and increasing correlations among particles that we expect with collapse and thus 
provide valuable insight into the nature of this instability. 

Let us now turn our attention to some analytical properties of the MFA. We will for the moment work with unitless 
parameters such that (3* — f3e, p* = pa d , and k* = ka. For notational convenience, we will continue to denote these 
quantities as /3, p, and k, respectively. It is important to note that when we speak of the infinite-dimensional limit, 
we mean the limit as d — > +00 such that [3* and p* arc held constant. Define the dimensionally-dependent parameter 
Xd = (3pir d / 2 , which implies (see (Eol) and (f42"|): 



S'MFA(fc) 



1 + A d exp[-^] 



(43) 



It is clear from (43) that Smfa -> 1 as i -> +00 for A^ fixed; indeed, one can see from Figure [5] that Smfa represents 
a "smoothed step function" for any finite dimension. 

Consider the limit Ad — > +00. Equation (43) implies that .Smfa — ► on any compact subset of M + (used here 
to denote the nonnegative reals) in this case. We note, however, that this limit does not commute with the limit 
k — ► +00; this issue is directly related to the apparent disappearance of the "step function" in the high-dimensional 
limit. To address this problem, it is worthwhile to consider the evolution of the slope of Smfa at Smfa(&) = 1/2, 




the "midpoint" of the structure factor, as — > +00. Assuming (3p > and fixed d, SmfaCO = 1/2 occurs at 
k c i = 2\J\n{\d) > 0. The slope of S'mfa at any point k in its domain is given by the derivative Smfa- 



Evaluating (45 1 at k^- 



>$mfa(^) 



Xdkexp 



2(1 + A d exp[~f]) : 



[S(k)Y 



fXdh\ 


" k 2 ' 




4 



•Smfa(^) 



k=k d 



(44) 
(45) 

(46) 
(47) 



We immediately notice from (47 1 that the slope diverges in the infinite-dimensional limit, indicative of the behavior of 
a true Heaviside step function. Our claim based on this information is that the MFA approaches a step function with 
discontinuity at k = +00 in the infinite-dimensional limit. We draw the connection here with a system of identical 
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FIG. 4: Hypernetted chain approximation to 172 with d — 5 demonstrating the (non-divergent) collapse phenomenon for value 
of /3e < 0. The density per 5 = 0.1. 



hard spheres, the pair correlation function of which is exactly a Heaviside step function; as d — -> +00, the mean field 
approximation to the GCM thereby resembles a system of hard spheres with arbitrarily large radii interacting in the 
dual space (in the sense of Fourier transforms) to the real space of the Gaussian core particles for nonzero density. 
In low-dimensional reciprocal space, the hard spheres are "smoothed" by the MFA, meaning interparticle penetration 
becomes increasing likely with decreasing dimension; however, the particles adopt an increasingly hard core as the 
dimension increases. This information allows us to draw some analytical conclusions about a decorrelation principle 
for the GCM in this regime. Assuming for sufficiently high dimension we may write as an approximation for S"mfa : 



<Smfa = 0(k-k d ), 



(48) 



where 0(fc — kd) denotes the Heaviside step function in reciprocal space with discontinuity at k = kd, it is possible to 
evaluate the coordinate space correlation functions analytically. The result is: 



dk 



g 2 {r) = l-v 



2tt 



d/2 



Jd/2(k d r) 



d/2 



(49) 
(50) 



where v — 1/ p denotes the reciprocal density. Since we consider d > 1, we utilize the principal asymptotic form of 
the Bessel function in (50 1 to obtain the result: 



92 (r) 



2 



d/2 



■Kd) 



(d+l)/2 



(51) 
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We make note of the constraint that (3p > n d l 2 to ensure that € K, restricting this approximation essentially 



to the low-temperature/high-density regime. Equation (51) is single- valued for all r; it is clearly less than 1 for any 
finite dimension, and the value it adopts is roughly the minimum of the expression for 92 in (50 1. Note that for 
d — > +00 the values of <?2 in (50) and (51) both approach #2(7) = 1. The fact that ( 50 1 and (51 1 converge in the 
limit d — ► +00 thereby reflects a loss of correlations in the infinite-dimensional limit, and from our knowledge of the 
behavior of the structure factor in the MFA, we expect that this "hard sphere" approximation accurately captures 
the high-dimensional behavior of the system. These results thereby provide analytical support for a decorrelation 
principle with the GCM in the fluid phase. 

To examine the behavior of the MFA in the high-temperature limit, we derive the cluster expansion of g^ 1 ^^ with 



respect to reciprocal temperature /3. 
respect to h. The result is: 



This is easily accomplished via iteration of the Ornstein-Zernike equation with 



(MFA) 
#2 



+00 



(r) = 1 - ^(r) + > (-(3) n C n+1 (r) 



(52) 



The MFA therefore keeps only the dominating terms from the analytical high-temperature expansion. However, 
unlike the analytical expansion, the series in (52 1 is convergent; the proof is left to Appendix D. Plots of gi obtained 
from (52 1 are given in Figure 6] For values of /? > such that the series converges, the form of <?2 is very similar 
to the results obtained from the cluster expansion work above. In fact, we expect that the MFA becomes a better 
approximation to the cluster expansion as d increases due to the domination of the series by the chain diagrams. 
Since the density expansions for the HNC and PY approximations, which are well-known (see, e.g., McQuarrie^) and 
are not derived here, retain the chain diagrams that appear in (52 1, we conclude that these approximation schemes 
should converge in the high-dimensional limit. For values of (3 < and within the radius of convergence, the MFA 
expansion shows a similar aggregation phenomenon to the one observed from the HNC approximation; as before, this 
approximation scheme penetrates into the instability region without collapse up to some finite value of /3, beyond which 
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FIG. 6: Radial distribution function from the series representation of the MFA in ( 52 1 for d = 5, 7 and (3 = 0.18, —0.18. Note 
that as the MFA passes through = 0, the form of g2 reflects the aggregation of the HNC approximation near the collapse 



instability; however, the series in (52 l is still convergent for some values of /3 < 0. 



the series diverges. We remark that the expression in (52 1 provides a computationally convenient way of computing 



approximations to 52 for the GCM to any order (3 n since analytic results are available for the chain diagrams. 



IV. SOLID-PHASE BEHAVIOR OF THE GCM 



A. High-dimensional lattice structures 



Our focus in the study of the solid phase of the GCM model will primarily involve the determination of lattice 
summation energies for known lattices up to d — 8. Unfortunately, structural information for lattices in dimensions 
significantly higher than eight is either unavailable or computationally prohibitive to obtain. However, the results 
presented here provide significant evidence in favor of a decorrelation principle. 

A lattice in high dimensions is characterized by the integer- valued linear combinations of a set of primitive basis 
vectors; i.e., a lattice A = {p}, where: 



E 



(53) 



In the expression above n, € M Vi, and denotes the i-th basis vector for the lattice. Associated with a lattice A is 
the so-called dual lattice A* with basis vectors q defined by (q, p) = 27rm, rn E Z. The lattice summation energy is 
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the total energy per particle for a given lattice, defined mathematically by: 



$(xi,...,xjv) = ^exp \-iCjaf] , (54) 



J=2 



where the a is the nearest-neighbor distance within the lattice and £j is a scaling factor that identifies the position 



of particle j relative to particle 1. Passing to the thermodynamic limit and partitioning the summation in (54 1 such 
that it is over all coordination shells in the lattice gives the desired result: 



The factor Z v denotes the coordination number of the i/-th coordination shell in the lattice. We note that the nearest- 
neighbor distance to the power d is inversely related to the density (i.e., a d p — c\), and the constant of proportionality 
depends on the chosen lattice. Therefore, all that is needed to completely specify the lattice summation energy for 
a given lattice are the sets {Z v } and along with the proportionality constant c\. We consider here the integer 
lattices Z d and the lattice families Dd,Ad, and Ed along with their respective duals. The Dd lattices are the d- 
dimensional counterparts to the three-dimensional FCC lattice and are the densest known lattice packings for all 
d < 6; similarly, the Ad lattice family generalizes the d = 2 triangular lattice. The Ed family contains the densest 
known lattice packings for d = 6, 7, 8 but are not uniquely defined for d < 6. 

B. The fluid-solid phase transition in the high-density limit 

From our knowledge of the three-dimensional phase diagram of the GCM, we expect to find a fluid-solid phase 
transition in any dimension d such that T m — > in the limits p — > and p — > +oo. The former limit is expected by the 
reduction of the GCM to a system of hard spheres in this regime, which has been shown with generality by Stillinger.^ 
There is strong numerical support (though no rigorous proof) for a fluid-solid phase transition with hard spheres up 
to d = 3, and we strongly suspect this is still true for higher dimensions. For d = 3, the freezing temperature for a 
hard-sphere system scales a^Tf(p) ~ exp (— cp~ 2 / 3 ), where c is a constant, and we conjecture that for arbitrary d 
the scaling is similar. In any case, our focus here is on the limit p — ► +oo. We initially consider a finite system of N 
particles in a volume f2 and introduce so-called collective coordinates, defined such that: 

$(x 1) ... J x J¥ )= # r «) ( 56 ) 

l<i<j<N 

l)^ n (k)[p(kM-k)-JV], (57) 



k 



where: 



6 n (k)= f 0(||x||)exp[-i(k,x)]dx (58) 

N 

p(k) = ^expH(k,x i )]. (59) 



i=l 



In passing to the thermodynamic limit, we have that <f>n — > 0, which is given in Q. 
Let ko be defined as a reciprocal space radius such that (3cf>(ko) = |; therefore, 



k Q = 2Jln(27r d / 2 (3). (60) 
The results in ^ imply that <f>(k) decreases as k increases. As a result, in the low-temperature regime for ||k|| < ko, 



(60 1 implies that f3<j)(k) is large and positive. By (57), minimization of the energy will then require p(k)p(— k) to move 
toward its minimum to offset the effect of increasing <p(k). However, for |jk|| > fc P4>(k) will become increasingly 
small, and the magnitude of p(k)p(— k) is of less consequence in the minimization of <f>. 
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Since the density of k- vectors inside a sphere of volume f2 is f2/(27r) d , the number Nq of such vectors with magnitude 
k < ko is: 



TT d / 2 k$ 



\T(d/2 + l)y \{2ir) d 7 2 d T(d/2 + l)7T d / 2 



(61) 



By the argument above, we conclude that for sufficiently high density p, Nq reflects the number of "lost" degrees of 
freedom to the system by the argument above. When Nq reaches some characteristic fraction < 8 < 1 of the total 
number Nd of degrees of freedom, the GCM will freeze; i.e.: 



(In [2(3 f ^ 2 ] ) d/2 n 



6Nd 



n d / 2 T[d/2 + 1] ' 

where /?/ is the reciprocal freezing temperature. Solving (|62| for ksTj yields: 



k B T f = 27r d/2 exp 



(d6fm d/2 T[d/2 + 1] 



2/d' 



(62) 



(63) 



which approaches in the limit p — > +oo. 



C. Duality relationships 

Since the Gaussian is self-similar under Fourier transform, it is possible to relate the lattice summation energy of a 
lattice at low density to the lattice summation energy of its dual lattice at high density; we call such an expression a 
type of duality relation. A dimensionally-dependent duality relation derived in this study for the lattice summation 
energies in the GCM is given below; however, we also mention the duality relationships recently put forth by Torquato 
and Stillinger regarding the ground state of a classical system interacting via a bounded, absolutely integrable pair 
potential </)(r)p2l One result is (based on the FT convention used in ([3])): 



,Ar)h(r)dr = [ — ) / < >[k)Ji[k )dk. 164) 



where h denotes the total correlation function as defined in (15 1, and h, (f> denote the Fourier transforms of h and (f>, 
respectively. Equation (64 1 is an immediate consequence of Parseval's formula (for a reference, see Lieb/Loss^J) since 



under the given assumptions (f>,h G L (R ); furthermore, one may show that if the configuration of particles in 



is a ground state and ergodicity is assumed, then the left- and right-hand sides of (64 1 are minimized. It should be 



stressed, however, that (64 1 will hold regardless of whether the configuration is a ground state. Equation (64 1 may in 
turn be used to prove the following duality relationship for a Bravais lattice Ap2l 

0(r = O)+ 0(r) = /#: = Q)+p ( 65 ) 

rGA\{0} keA*\{0} 

which follows from the identity h(r) = X) Z n 8{r — r n ) — 1 for a Bravais lattice, where si(r) denotes the surface 

n— 1 

area of a d-dimensional sphere. 

Torquato and Stillinger go on to show that twice the minimized energy per particle for any ground-state 

structure of the dual potential <f>(k) is bounded from above by the corresponding real-space minimized twice-energy 



per particle <I> m i n , i.e., the right-hand side of (65 1 



p4>(k = 0)+p ]T $(k). (66) 

k£A*\{0} 

This inequality results from the notion that the energy-minimizing configuration in the dual space to a real-space 
configuration need not be a Bravais lattice. At the very least, such a possibility cannot be eliminated solely from 



(651. However, equality of minimum energies in real and reciprocal spaces will hold whenever the reciprocal lattice 
A* at reciprocal lattice density p = p~ 1 (2n)~ d is a ground state of <j){k). Conversely, if a sufficiently well-behaved 
dual potential 4>{k) has a Bravais lattice A* at number density p, then: 

p(b(r = 0) + p ]T <P(r). (67) 

r£A\{0} 
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Here we present a duality relationship that associates the energy per particle ($/N) A of a given lattice A in the 
thermodynamic limit at low density with the equivalent energy per particle of the dual lattice A* at high density. 
In accordance with prior work by Stillinger^Sl concerning one, two, and three dimensional duality relations for the 
GCM, we consider the energy per particle ($/iV) to eliminate boundary effects in passing to the thermodynamic limit. 
Define: 



/2<£ \ 

I A (a) = 1 + lim — , 



(68) 



where a denotes the nearest-neighbor distance within the lattice. We may equivalently write I\(a) in terms of the 
discrete density function g(s): 



I A (a) 



/ g(s) exp(— s 2 )ds, 
Js. d 



where g(s) is given in terms of a summation over Dirac delta functions: 

0(S) = £<^>(S- Sj ). 



(69) 



(70) 



It is convenient at this point to "smooth" the Dirac delta functions in (70 1 via convolution with a normalized Gaussian, 
yielding: 



g(s) = lim g(s, a) 

a — >+oo 

£>(s,a)=(-) 1 J^exp [-a(s - Sj) 2 ] , 



(71) 
(72) 



where it is understood that the limit a — > +oo is to be taken at an appropriate point in the calculation. Noting that 
g(s,a) is a d th -order periodic function of the variable s, we may represent this function in a Fourier series: 



Q{s, ot) = ^2 /( k ) ex P b'( k > s )] 



(73) 



where (x,y) = Xjfji denotes the inner product of two vectors in d-dimensional real Euclidean space, and the 

vectors k are 2tt times the vectors from the dual lattice A*. Expressions for each /(k) are determined from Fourier 
orthogonality conditions, whereby one multiplies ( 73 ) by exp [— i(k, s)] and integrates over a unit cell within the lattice 
to obtain: 



/(k) = p A (o) exp 



-k 



21 



4a 



(74) 



where pa (a) = c^/a d denotes the density of the system as a function of the nearest-neighbor distance a with ca a 
constant for the lattice A. Combining (71 l-(74l yields: 



Ia(cl) — lim / /(k) exp [i(k, s)l exp(— s 2 )ds 
= lim > /(k) / exp [i(k, s)] exp(— s 2 )ds 

a^+co J Rd 



n d / 2 pA(a) lim \^exp 

rv- — >4-rxn ' ^ 



-fc 2 (q + l) 
4a 



TT d/2 PA (a)J2 e ^P 



-k 2 



(75) 
(76) 
(77) 
(78) 



The right-hand side of (78 1 is exactly of the form for 1(a) given in (68 1; in fact, under suitable scaling for the dual 
lattice A*, we may write: 



I A (a) = 7r d / 2 p A (a)5>xp{- {ir X ub(a)} 2 } 



(79) 
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where b(a) denotes the nearest-neighbor distance for the dual lattice as a function of a, and %k denotes the related 
scaling factor for the particle coordinate in the dual lattice. Reference to (68 1 yields the desired result: 



I A (a) = n d / 2 p A (a) J A .[7r6(o)]. 



(80) 



Our interest here is to discern the so-called self-dual density p* , which is the density at which the lattice A and its 
dual A* have the same energy per particle in the thermodynamic limit. It is immediately clear from (80) that if the 
condition of equal energy per particle between a lattice and its dual holds, then: 



P 



-d/2 



(81) 



whereby the coefficient on the right-hand side of (80 1 becomes unity. 



D. Lattice energies and coexistence regions 



The lattice summation energies for our chosen lattice families (Dd, Ad, Ed) relative to the energy for the corre- 
sponding 7L d lattice are given in Figures In accordance with our predictions from the duality relationship 
in (80), there appears to be a universal phase-transition density between and lattice and its dual at p — -K~ d / 2 . For 
d = 4, 5, the lowest-energy lattices are given by the Dd lattice and its dual; this observation provides direct support 
for the Torquato-Stillinger conjecture^ concerning the ground-state structures of classical systems since Dd is the 
Bravais lattice corresponding to the densest known sphere packing for 3 < d < 5. As we would predict from this 
conjecture, the Ed lattices and their duals obtain the lowest energy for 6 < d < 8, followed by the Dd family. Thus, 
there appears to be a relationship between the density of a given lattice structure and its energy with respect to the 
GCM, providing numerical support for the conjecture mentioned above. 

We may define the lattice coexistence region between dual lattices via the Maxwell double-tangent construction. 
For notational convenience, we will define: 



e = ($/AT) 

v = (1/p) = (V/N). 



For a one-component system, the internal energy $ satisfies the following differential: 

d$ = TdS - pdV + pdN, 



(82) 
(83) 



(84) 



where p is the chemical potential of the system. Since we will be interested here in the variation of the internal energy 
with reciprocal density, we will without loss of generality assume constant entropy and simply write: 



d$ = ~ p dV + pdN. 



(85) 



The Maxwell double-tangent construction involves finding a solution (v^- A \v^ A )) to the following set of coupled 
equations: 



p(A)(«(A))=p( A *)(D( A *)) 

M (A)( S (A) )=M (A*) (s (A-) ); 



(86) 
(87) 



where and v^ A ) characterize the upper and lower reciprocal density bounds for the phase coexistence region 
between a lattice and its dual. That ( 86 1 and ( 87 1 correspond to a double-tangent is an immediate consequence of 
Euler's theorem for first-order homogeneous functions (a succinct review may be found in Chandler- 8 ) and (85 1, which 
imply: 



-pV + pN. 



Therefore, (87 1 is equivalent to: 



(e + pv) ^ = (e + pv) 



(A*) 



(89) 



The equal-pressure condition in (86 1 allows us to write the following equation for a line tangent to e' a ' and e' A > at 
and y( A \ respectively: 



r (A) _ „(A*) 



-,(A) 



-,(A*) 



(90) 
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FIG. 7: Relative lattice summation energies (&/N) zd — (&/N) A v. 1/p for given lattice families A^; d — 4,5,6. 



To solve for and v^ A \ we utilize (85) to obtain 



P = 



V = 



A'* 



9$ 
dV 

dNj v ~ { dN ) 



= (d(Ne)\ 



N 



de 
ON 



N 



vp. 



dv 
dN 



£ — V 



(91) 

(92) 

(93) 
(94) 



Equations ( |91| and (94) allow us to numerically determine values for u^ A ' and v^ A > from our lattice summation data. 

Table |l] collects the results for the Maxwell construction of the phase-coexistence regions between dual lattices in 
each dimension. We note that the width of the phase coexistence region scaled by the self-dual density p* = 7r~ rf / 2 
increases with dimension for a given lattice family. The immediate significance of this behavior is that for arbitrarily 
high Euclidean dimensions, it is possible that the phase coexistence region for a particular lattice is wide enough such 
that new structures are able to achieve lower energy as ground states. In the context of a decorrelation principle, we 
cannot exclude the possibility that these structures are disordered. If the conjecture by Torquato and Stillingep^l that 
the densest known packings of hard spheres in high dimensions are disordered is to be believed, then this behavior is 
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FIG. 8: Relative lattice summation energies ($/JV) zd — (3>/iV) A v. 1/p for given lattice families A<j; d = 7,8. 
expected. 



d 


Lattice Family 


P* 


P* 


PA 


PA* 


\pA - PA*\/p* 


4 


A 


0.47983 


0.03310 


0.10103 


0.10164 


0.00605 


5 


A 


0.53922 


0.01980 


0.05683 


0.05749 


0.01154 


5 


D 


0.51155 


0.01927 


0.05707 


0.05726 


0.00332 


6 


A 


0.60331 


0.01186 


0.03196 


0.03254 


0.01783 


6 


D 


0.56598 


0.01146 


0.03206 


0.03244 


0.01154 


6 


E 


0.54849 


0.01136 


0.03236 


0.03244 


0.00248 


7 


D 


0.62835 


0.00684 


0.01800 


0.01841 


0.02257 


7 


E 


0.57723 


0.01826 


0.01820 


0.01826 


0.00330 


8 


D 


0.69825 


0.00410 


0.01009 


0.01045 


0.03483 



TABLE I: Scaled phase coexistence regions for specified lattices in 4, 5, 6, and 7 dimensions. Here p* and p* are the self-dual 
chemical potential and pressure, respectively; p* is the self-dual density as defined in (pTJ above. 



V. CONCLUDING REMARKS 



We have hereby made an effort to generalize the phase properties of the GCM up to d = 3 with respect to the 
dimensionality of the system. In the fluid phase, we have developed a dimensionally-dcpendent low-density /high- 
temperature expansion for the radial distribution function gi- Although this series provides some interesting evidence 
for a decorrelation principle, it suffers from the drawback of being divergent due to a collapse instability induced by 
the unmediated attraction of the particles for values of (3 < 0. Numerical approximations (HNC, PY, and MFA) for 
gi at higher densities and lower temperatures strengthen the evidence for decorrelation in the GCM, and we show 
that the MFA contains the high-dimensional behavior of the GCM for sufficiently high temperatures. With regard 
to the solid phase of the GCM, results for the lattice summation energies of known lattice families up to d = 8 
provide direct support for the recent Torquato-Stillinger conjecture^ concerning the ground states of classical many- 
particle systems. Namely, for sufficiently well-behaved pair potentials the ground states correspond to the densest 
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known lattice packings for low densities and the corresponding dual lattices for high densities with a solid-solid phase 
transition at intermediate density values. Having mentioned this idea, we cannot exclude the possibility that for 
sufficiently high dimensions the ground-state structures may in fact be disordered; in the case where the particles 
adhere to a decorrelation principle, it is not unreasonable that we might expect this result in accordance with prior 
work by Torquato and Stillingcr. 9 

Despite our current work on this problem, a few remaining points are worthy of mention. With regard to the 
behavior of the system in the fluid regime, the nature of the collapse instability makes analytical evaluation of 
the radial distribution function difficult to interpret. We have developed the density expansion of gi and made 
favorable comparisons to numerical approximations, yet the physicality of the cluster expansion seems to be lost due 
to divergence of the series. Stillingei 17 priorly utilized a Borel resummation to interpret the "lost" information in 
the high-temperature expansion of the excess free energy, yet the strong spatial and dimensional dependence of the 
gi expansion makes this technique difficult in the present case. What is perhaps most promising is the reduction 
of the expansion to that of the MFA in the high-dimensional limit, whereby the divergence of the series is partially 
removed with results similar to the numerical approximations from the HNC. However, the MFA has been shown 
to be a thermodynamically inconsistent approximation^ meaning that it cannot be an exact theory to describe the 
GCM. 

In the solid phase, we still have an open question concerning the fluid-solid phase transition at any value of p. There 
is reason to believe, as mentioned above, that such a transition exists up to some maximum melting temperature T m (p) 
as in the three-dimensional caseP and we have implicitly made this assumption in the present work. However, we 
know that the phase diagram will be dimensionally-dependent as shown by the fact that for d = 4, the self-duality of 
D4 preempts a solid-solid phase transition as in d = 3. Nevertheless, the study of the lattice structures of the GCM 
is worthwhile in the context of the Torquato-Stillinger conjecture^ concerning the ground states of classical systems 
and with respect to the decorrelation principle. Though the proof of the existence of disordered ground states for the 
GCM in high dimensions is still unavailable, our results show that we certainly cannot preclude the possibility. 

The results we have presented suggest that it may be possible to extend this research to a broader class of pair 
interactions. We recall that the primary advantage of the GCM, as we have shown, is the property of being self- 
similar under Fourier transform. However, it is known that the cigcnfunctions of the Fourier transform are equivalent 
to the states of the quantum- mechanical harmonic oscillator (under suitable scaling), namely Gaussian functions 
modulated by the Hcrmitc polynomials. This property indicates that our analysis may be generalized to pair potentials 
containing both repulsive and attractive components. This more general case is likely to modify the collapse instability, 
and determining the precise manner of this modification is an area that deserves further study. Physically, these 
interactions could be related to spatially inhomogeneous solvent compositions that simultaneously induce repulsion 
and attraction among macromolecules in solution. Alternatively, one may also consider pair potentials formed from 
linear combinations of (attractive and repulsive) Gaussian interactions. In any case, the possibility for finding unique 
thermodynamic phenomena which could be used in the design of novel materials makes this a viable avenue for future 
exploration. 
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APPENDIX A: PROOF OF REGULARITY OF THE GCM PAIR POTENTIAL 
Proposition: The GCM pair potential is regular for e > 0. 

Proof: From the definition of regularity, it is sufficient to prove the proposition for (3 > 0. Fix /3, e > 0. Positivity 
of the pair potential Q provides the requisite lower bound of K = B = 0, where B denotes the stability constant for 
the pair interaction. Since < exp [— /3</>(x)] < 1 for /3, e > 0, we thus have: 

C(fl)= [ {l-exp[-/rexp(-ax 2 )]}dx, (Al) 



where a = 1 /a 2 and (3* — (3e > 0. 

Since the exponential function is an entire function, we may write the integrand of ( Al I as a uniformly convergent 
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Taylor series: 



(A2) 



The series in (A2| is uniformly convergent on any compact interval in R, which implies we may exchange the 



operations of integration and summation to obtain: 
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(A3) 
(A4) 



One may show that the series in (A4) is convergent via reference to the ratio test. 

d/2 



lim ( n = lim 

n— *+oo n— H-oc 



= /3* lim 



(/3*) n 



1 



(n+1)! 



71+1 



n— >+oo I \n + 1 J \n + 1 
= < 1 



(A5) 

(A6) 
(A7) 



Convergence of the series in (A4) proves the proposition. □ 



APPENDIX B: PROOF OF EQUATION (32) 



Proposition: The chain diagrams defined by (31) satisfy: 



C n +i — p" 



-IN d/2 



(Bl) 



Vrt > 2, where denotes the GCM pair potential. 

Proof: We prove by mathematical induction. Let p(k), k E N\{1}, be the claim: 



/ 7T 

U 0*,(i+i)dx 2 ' ' ' rfx ^ = ( ~~ 



r k-l\d/2 



(B2) 



As mentioned in the text, we have chosen the stationary vertices in the cluster diagram to be 1 and k+ 1 for simplicity. 
To prove p(2), we use the relation: 



m + n 



d/2 



m + n 

.k 1 



(B3) 



which follows directly from ([8]) 
Equation (B3) implies: 



01,202,3^X2 = (I) <t>l,t, 



(B4) 



which is the exact form of the right-hand side of (Bl I. Therefore, p(2) is true. 
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Fix ko £ N \ {1} arbitrarily, and make the usual induction hypothesis that p(fco) is true. We now prove p(ko + 1) 
is true. Using Fubini's theorem: 



/feo+i p 
Yl ^i,(i+l) dx 2 ' ' ' dx ka + l = / i 
i=l ^ 



(fc +l),(fco+2) ( n^.(<+ 1 ) dX2 "" dXk ol - rfx fco+l 



(/c + l),(feo+2) ( ? !, ^(t", + l) dx feo + l 



x d/2 

\ ,l/(feo + D 



(B5) 

(B6) 
(B7) 



We see from the result in (B7l that p(ko + 1) is true, and the truth of p(k) Vfc G N \ {1} immediately follows. The 



proof of (Bl I is thus apparent with the n — 1 factors of p arising by convention from the integration over the n — 1 



integrable nodes. □ 



APPENDIX C: PICARD ITERATION ALGORITHM FOR NUMERICAL APPROXIMATIONS 



In order to obtain a numerical approximation to 52 using one of the closures mentioned in Section II, we are required 
to solve the OZ equation, which for a given closure c(r) will be a nonlinear integral equation. To address this problem, 
we first make an initial estimate for the convolution integral 7(7") = h(r) — c(r) and calculate the corresponding c(r) 

using either (35 1 or (36 1. To update our guess for the form of 7(7*), we use Y(r) = S^ 1 jp [c(fc)] 2 /[l — pc(fc)] j, 



where 



denotes the inverse Fourier transform, and the right-hand side of the equation follows from taking the FT of the 
OZ equation. We calculate the error e = \f&r\\~f' (r) — 7(r)||, where 6r denotes the size of the coordinate-space mesh, 
and set 7(7*) = 07' (r) + 0- ~ a )"fi r )-> where < a < 1 is a mixing parameter to speed (or to aid) convergence of 
the algorithm. The algorithm is iterated so long as e is greater than some specified tolerance. Convergence of this 
algorithm is reasonably quick assuming a good choice of a, requiring usually only a few hundred interations. We 
note, however, that as the density increases, it is necessary to decrease the value of a, thereby increasing the requisite 
iterations. 



APPENDIX D: PROOF OF CONVERGENCE OF THE SERIES IN (52) 
Proposition: The function g^ FAS> is well-defined by (52 1; more specifically, 



+00 



^(-/?)"C n+1 (r)<+oc 



TJ=2 



Vr G [0, +00) within some finite radius of convergence with respect to (3. 

Proof: We show that the series in (Dl I converges by reference to the ratio test. Namely, 
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(D2) 

(D3) 
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Therefore, the series in ( |DT) will converge for all /3 such that \/3\ < l/(pir d / 2 ). (Note that we have, without loss of 



generality, chosen unitless parameters as described in Section III.) □ 
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